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We present an exact numerical stndy of the scaling of density and momentum distribution func¬ 
tions of harmonically trapped one-dimensional bosons with repulsive contact interactions at zero 
and finite temperatures. We use path integral quantum Monte Carlo with worm updates in our 
calculations at finite interaction strengths, and the Bose-Fermi mapping in the Tonks-Girardeau 
regime. We discuss the homogeneous case and, within the local density approximation, use it to 
motivate the scaling in the presence of a harmonic trap. For the momentum distribution function, 
we pay special attention to the high momentum tails and their k~'^ asymptotic behavior. 
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I. INTRODUCTION 

Developments in ultracold atomic gas experiments 
have boosted the study of many-body interaction effects 
in bosonic systems m- In particular, the ability to load 
Bose-Einstein condensates in optical lattices has provided 
a unique opportunity to control the effective dimension¬ 
ality and interactions of Bose gases. A setup that is of 
interest to us here is that of a BEC loaded in a two- 
dimensional (2D) optical lattice, as a result of which it 
splits into an array of cigar-shaped Bose gases at very 
low temperatures [5]. For a sufficiently deep 2D opti¬ 
cal lattice and low energies, each cigar-shaped Bose gas 
can behave as an effective one-dimensional (ID) system 

mm that is described by the Lieb-Liniger model mm 

(plus an additional term to account for trapping poten¬ 
tials). The effective ID interaction strength is deter¬ 
mined by the three-dimensional (3D) scattering length 
and the transverse confinement provided by the lattice 
[7]. For sufficiently strong effective ID interactions and 
sufficiently low ID densities, these gases enter the so- 
called Tonks-Girardeau regime Eli, in which bosons 
behave as impenetrable particles (hard-core bosons) and 
can be mapped onto noninteracting fermions |10j . Fur¬ 
thermore, if an additional (weaker, but not too weak) lat¬ 
tice is added along the longitudinal direction, the system 
can be described by the ID version of the Bose-Hubbard 
model, which has been of much interest in condensed 
matter physics m- Experiments in the presence of that 
additional lattice have allowed for the observation of the 
superfluid to Mott insulator transition in one dimension 
[T^ and the lattice version of the Tonks-Girardeau regime 

m- 

The Lieb-Liniger model has attracted quite some at¬ 
tention since its introduction more than 50 years ago 
Hill]. It is an integrable model whose exact solution, 
which can be obtained using the Bethe ansatz m , pro¬ 
vides insights into the universal behavior of ID gapless 
systems. Its far-from-equilibrium dynamics has been re¬ 
cently scrutinized theoretically [T5IBD] , after experiments 
showed that ID Bose gases taken far from equilibrium 
relax to states in which observables are not described 
by traditional ensembles of statistical mechanics m and 


generalizations of the latter are needed [HI |H] ■ 

We should stress that despite the fact that the Lieb- 
Liniger model is Bethe ansatz solvable, the calculation of 
correlation functions is extremely challenging due to the 
complexity of the Bethe ansatz eigenfunctions. In the 
ground state, for example, the evaluation of one-particle 
correlations (of which the momentum distribution func¬ 
tion is the Fourier transform) has only been possible nu¬ 
merically within algebraic Bethe ansatz [H] , and through 
diffusion Monte Garlo simulations [HI ES] ■ At finite tem¬ 
peratures, they have been calculated using simulations 
within the stochastic gauge method [27|- Remarkably, 
the asymptotic behavior of the momentum distribution 
function [m(fc)], for large values of the momentum k (as¬ 
sociated with short range correlations), has been deter¬ 
mined analytically [m{k) oc k~‘^] [28] ■ In the Tonks- 
Girardeau regime, the simplification introduced by the 
mapping to noninteracting fermions has enabled the an¬ 
alytic calculation of one-particle correlations at zero and 
finite temperatures [HHH] . On the other hand, the ID 
Bose-Hubbard model is not integrable, so correlations 
can only be computed numerically [2|- However, its 
Tonks-Girardeau limit is integrable and can be mapped 
onto the so-called XX chain [5], which is the isotropic 
limit of the XY chain introduced by Lieb, Shultz, and 
Mattis [33j . and which can be mapped onto noninteract¬ 
ing fermions [2] . This way, it is possible to calculate ana¬ 
lytically one-particle correlations of the Tonks-Girardeau 
gas in the presence of a lattice EH ES] • 

In experiments with ultracold gases, a confining po¬ 
tential (that is, to a good approximation, harmonic) is 
needed to contain the gas. In the presence of such a 
confining potential, as in the Hubbard model, only the 
limit in which bosons are impenetrable (hence, mappable 
to noninteracting fermions) remains integrable. In this 
limit, the ground-state density and momentum distribu¬ 
tion functions, as well as one-particle correlations, have 
been studied in the continuum |3tiH41] and in the lattice 
[42| [43]. In those studies, universal power laws with the 
same exponents as the ones known from homogeneous 
systems were found. At finite temperatures, a system¬ 
atic study of one-particle observables has only been re¬ 
ported in the lattice [44] . while, in the continuum, the 
momentum distribution function was recently calculated 
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in Ref. |45| for trapped systems with five particles. For 
finite interaction strengths, for which trapped systems 
are not integrable anymore, quantum Monte Carlo sim¬ 
ulations (with a lattice discretization) have been used to 
compute momentum distribution functions in the con¬ 
tinuum and compare them with experimental results at 
relatively weak interactions, obtaining a good agreement 
between the two | 1 S] . 

In this work, we are interested in systematically study¬ 
ing the density and momentum distribution functions of 
the Lieb-Liniger model in a harmonic trap at zero and 
finite temperature. We will focus on the scaling proper¬ 
ties of those quantities and in the high momentum tails 
of the momentum distribution function. The latter have 
remained elusive to quantum Monte Carlo simulations 
so far. The weight of the high momentum tail, which is 
known as Tan’s contact, is a quantity that plays a central 
role in a set of universal thermodynamic relations known 
as Tan’s relations [T7ll49| . The Tan’s contact has been 
measured in experiments with 3D Bose-Einstein conden¬ 
sates m- 

In order to obtain very accurate results for all the quan¬ 
tities above, for arbitrary interaction strengths, here we 
use path integral quantum Monte Carlo |51j for contin¬ 
uous systems with worm updates [SH . We also ob¬ 
tain results in the Tonks-Girardeau limit using the Bose- 
Fermi mapping in the lattice and working at very low 
fillings For completeness, in all cases, we also 

discuss the behavior of the quantities of interest in ho¬ 
mogeneous systems, which helps motivate the scaling re¬ 
lations for the trapped systems. 

The paper is organized as follows. In Sec. [TJ we intro¬ 
duce the Hamiltonian and numerical approaches used. 
We also discuss some of the checks done to gauge the 
accuracy of our calculations. In Secs. [Tn] and |IV) we 
report a detailed study of density and momentum dis¬ 
tribution functions, as well as one-particle correlations, 
of Lieb-Liniger systems at zero and finite temperatures, 
respectively. In both sections, we discuss results for ho¬ 
mogeneous and harmonically trapped systems. Finally, 
our conclusions are presented in Sec. |V] 


II. HAMILTONIAN AND NUMERICAL 
APPROACHES 


We consider ID bosons with repulsive contact interac¬ 
tions in the presence of an external harmonic trap. The 
Hamiltonian can be written as 
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2m dx'j 


+ V{xj) 


+ 9'^Hxj - xi), 
j<^ 


( 1 ) 


where m is the mass of the bosons, g is the strength 
of the contact interaction, V{xj) = ma;^a;|/2, and w 
is the frequency of the harmonic trap. The strength of 
the effective ID contact interaction is usually written as 
g = —2h^/maiY^, where uid = —a_L(a_L/as — C) is the 



FIG. 1. (Color online) Worm algorithm results for the ground- 
state kinetic and interaction energies (see text) of the Lieb- 
Liniger model (circles) compared with the exact analytical 
results (lines) [3. The statistical errors (shown within the 
circles) can be seen to be much smaller than the circle sizes. 
The average number of particles in the system is Nb ~ 20. 


ID scattering length, Os is the (3D) s-wave scattering 
length, a_L is the length of the transverse confinement, 
and C = 1.0326 [7]. Another parameter frequently used 
to describe the interaction strength is c = mg/h'^. In 
the absence of the trap, this Hamiltonian reduces to the 
Lieb-Liniger model, which essentially has only one tun¬ 
able parameter 7 = c/p, where p is the density. 

Here, in order to simulate the Hamiltonian above in the 
continuum, we use a path integral quantum Monte Carlo 
method with worm updates (from now on referred to as 
the worm algorithm) |521153j . The configuration space 
in our simulation is built with discrete imaginary time 
world lines in a continuous position space. The worm 
algorithm operates both with diagonal and off-diagonal 
configurations by introducing one additional open world 
line in the configuration space. This makes it very ef¬ 
ficient in simulating both diagonal observables (such as 
the energy, the density distribution, and the superfluid 
fraction) and off-diagonal observables (such as two-point 
one-particle correlation functions). In our calculations, 
the action for contact interactions is approximated by 
the pair-product action, see Appendix]^ and we always 
work in the grand canonical ensemble. 

In Fig. we compare results for the ground state 
kinetic [K{"f) = h?t{"f)/2m\ and interaction [U{j) = 
/2m\ energies from the worm algorithm simula¬ 
tions (for a finite system with « 20 bosons) with 
the exact analytic results in the thermodynamic limit [S] . 
Average energies are obtained using the thermodynamic 
estimators reported in Appendix (see also Ref. EH). 
Figure shows that, already for systems with iVf, « 20, 
finite-size effects are negligible, and that the worm algo¬ 
rithm provides a very accurate estimation of observables 
in diagonal configuration space. 

The observables on which we focus in this work, i.e., 
the density and momentum distribution functions, are 
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FIG. 2. (Color online) (a) Ground-state one-particle corre¬ 
lations and (b) momentum distribution functions of bosons 
for 7 = 8.33. Results are reported for systems with the same 
density but different average number of particles. Dashed 
lines are the results obtained using the worm algorithm and 
solid lines the results obtained in numerical calculations via 
algebraic Bethe ansatz |24) . 


calculated from the one-particle density matrix 

g{x,y) = {¥{x)'l’{y)), ( 2 ) 

where fE'^(a:) [^'(a:)] is the bosonic field creation (annihi¬ 
lation) operator at position x. g{x, y) is an example of an 
off-diagonal observable that can be efficiently calculated 
using the worm algorithm. 

For improved statistics, in our calculations, we average 
g{x,y) over all possible translations {0 < x,y < L) so 
that it becomes a function of r = |a; — y\. For r > 0, 


s(7 = j 


pL pL 

/ g{x,x + r)dx + / g{y + r,y)dy 
JO Jo 


, (3) 


while, for r = 0 , 



g(x,x)dx. (4) 


The Fourier transform of g(r) gives the momentum dis¬ 
tribution function 


m( 


i(fc) = ^ '‘h[x,y) 

L 

= — / cos{kr)g{r)dr. (5) 

27r Jq 

In Fig.[^ we show results for g{r) [Fig.j^a)] and m{k) 


[Fig.j^b)] for the ground state of systems with 7 = 8.33, 
the same density, and different average number of parti¬ 
cles, and compare them with numerical results obtained 
via algebraic Bethe ansatz [21]. The results for g{r) 
in both approaches are essentially indistinguishable in 
Fig. [^a). The same is true for m{k) in Fig. [^b), ex¬ 
cept at the highest momenta. In the latter regime, m(A:) 
becomes very small and the statistical errors in the quan¬ 
tum Monte Carlo simulation become of the same order. 
This is why it is a challenge to study the asymptotic 
behavior of m(k) at high momenta. Experiments will 
suffer from exactly the same limitation as m{k) will at 
some point become of the same order of the experimental 
noise. The results in Fig. make apparent the accuracy 
of the worm algorithm for computing off-diagonal one- 
particle observables, which are the focus of this work. 
Once translational invariance is broken, no analytical so¬ 
lution is available. To gauge the accuracy of the worm 
algorithm in that case we focus on the Tonks-Girardeau 
limit. 

In that limit, the one-particle density matrix (and, 
hence, the momentum distribution function) in a lattice 
can be calculated exactly via the Jordan-Wigner trans¬ 
formation and making use of properties of Slater determi¬ 
nants, both at zero [221133] and hnite temperature [ 33 ] . 
As we explain below, using this approach in the low- 
density limit in the lattice, one can efficiently compute 
one-particle properties in the continuum |43j . 

The Tonks-Girardeau Hamiltonian in the lattice can 
be written as 


where b] ( 6 j) is the creation (annihilation) operator of 
a hard-core boson at site i, rii = b\b^ is the site occu¬ 
pation operator, t is the hopping amplitude, and V sets 
the strength of the harmonic trap. This Hamiltonian is 
obtained from the Bose-Hubbard Hamiltonian by taking 
the limit in which the onsite repulsion [/ —>■ cx) [ 2 |. 

The hard-core boson creation and annihilation opera¬ 
tors satisfy standard bosonic commutation relations with 
the constraint {b\Y = _^ = 0. Mapping the Tonks- 
Girardeau Hamiltonian (m onto a spin-1/2 chain [5Hl55j . 
and then the spin- 1/2 chain onto noninteracting fermions 
in ID via the Jordan-Wigner transition, one gets 

= -tY.Ulh+1+^-c.) + (7) 

i i 


where // (/J is the creation (annihilation) operator of a 
spinless fermion at site i, and n{ = is the fermionic 
site occupation operator. The fermionic Hamiltonian 
([^ can be straightforwardly diagonalized. One-particle 
bosonic correlations can then be obtained using proper¬ 
ties of Slater determinants as discussed in Refs. [42ti44] . 
We should stress that, within this lattice approach, the 
ground-state calculations are done in the canonical en- 
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FIG. 3. (Color online) Density profiles (a) and momentum 
distribution functions (b) of the Tonks-Girardeau gas as cal¬ 
culated using the lattice approach (red solid lines) and in the 
continuum using the worm algorithm (dashed lines). The av¬ 
erage number of particles is 5. Results are presented for dif¬ 
ferent temperatures: T = ksT/huj=0.1, 5, 10. Thin dotted 
lines indicate the asymptotic behavior of the momentum 
tails. The statistical errors in the worm algorithm are of the 
order of the fluctuations seen in the results. They are not 
reported for clarity. 



FIG. 4. (Color online) Tan’s contact for the trapped sys¬ 
tems in Fig. 1^ compared with results in Ref. [1^. The 
contact from the worm algorithm simulation is obtained 
through the interaction energy using the expression Coho = 
—a%^QNi,v{'y)/{'RaiY)) [derived from Eq. (lOlj, where v{'y) 


is the average interaction energy per particle [in units of 
fi^/(2m)]. We calculated Ca^o for two small values of Uid 
obtaining results that agree with the ones in Ref. [45]. The 
contact from the lattice calculation is obtained via a linear 
fit of the high momentum tails, and then extrapolated to the 
zero density limit via a quadratic fit to obtain the result for 
zero a/uHo (see inset). 


In the lattice, the momentum distribution function 
is defined as the discrete Fourier transform of the one- 
particle density matrix 


semble [421143] , while finite-temperature ones are done in 
the grand-canonical one jUj. 

To establish the relation between the parameters in 
the lattice Hamiltonian and in the continuum, we notice 
that the single-particle energy spectrum of the spinless 
fermion model without a trap is = — 2 t cos(fca), where 
a is the lattice spacing. In the zero-density limit, the 
Fermi momentum kp —>■ 0. As a result, the energy spec¬ 
trum becomes quadratic in k as, up to a constant, we 
can write £). = ta^k'^. Thus, in the zero-density limit, 
the lattice model reduces to a continuum model with the 
effective mass given by m = k?/2ta^. In this limit, ac¬ 
cording to Eq. ([^, the external trapping potential is 

V{x) = Vix/af. (8) 

The trapping frequency is then given by w = 2\/fV //i, 
and so the harmonic oscillator length oho = 
can be written as ano/o = (t/F)^^^. Finally, in finite- 
temperature calculations in the lattice, the lattice tem¬ 
perature Tt is usually defined as tTt = kpT, i.e., it is 
given in units of t. The temperature in units of trap¬ 
ping frequency in the continuum is related to the lattice 
temperature by the expression ksT/Huj = Tt^Jt/V / 2 . 
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where N is the number of lattice sites. By comparing 
Eq. (§ and Eq. (|^, one can see that the momentum 
distribution function in the continuum and in the lattice 
are related by the expression ruk = m(k)5k, where 5k = 
2 pIL (with L = No) is the discretization of k in the 
lattice model. 

With the relations discussed so far at hand, we are 
ready to compare results for density and momentum 
distribution functions in the continuum and the lattice, 
which, as mentioned before, help us gauge the accuracy 
of the worm algorithm for diagonal and off-diagonal ob¬ 
servables in the absence of translational invariance. In 
order to approach the Tonks-Girardeau limit in the con¬ 
tinuum model, we choose a very small value of the ID 
scattering length Oid/oho = 5.66 x 10“^. 

The scaled density {n(x) = g(x,x)] and momentum 
distribution functions obtained in the lattice and the 
continuum are shown in Fig. at various temperatures 
(T = kBT/huj). The agreement between the results of 
both approaches is almost perfect at the lowest tempera¬ 
tures T = 0.1 and 5 (for T = 0.1 the system is essentially 
in the ground state). As the temperature increases, lat¬ 
tice effects become important and lead to visible differ- 
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ences when the lattice results are compared with the re¬ 
sults in the continuum [see the discussion after Eq. ®]. 
In the momentum distribution function [Fig. [^b)], the 
worm algorithm at T = 0.1 and 5 already captures part 
of the momentum tails. For higher temperatures, 
such as T = 10 in Fig. [^b), the k~^ regime is not visible 
in the worm algorithm calculation. This is understand¬ 
able because, as the lattice calculation clearly shows, the 
tail starts at higher values of the scaled momentum 
[smaller values of the scaled m{k)\ as the temperature 
increases. 

We have studied the weight C (Tan’s contact) of the 
k~'^ momentum tails [m{k) = C k~'^] using both methods. 
In the lattice approach, C is computed directly via linear 
fits of the momentum tails. In the worm algorithm cal¬ 
culation, the contact C is calculated via the interaction 
energy estimator, using the relation |45j 

( 10 ) 

The results are shown in Fig. The values of the con¬ 
tact obtained using the worm algorithm for two small 
values of the scattering length are consistent with the re¬ 
sults in Ref. [JS]. This indicates that such small values 
of oiD allow us to obtain results in the Tonks-Girardeau 
limit. Using the lattice approach, the contact in the con¬ 
tinuum is calculated by taking the low density limit by 
means of a quadratic fit (see the inset in Fig. |^. The 
results of such a fitting procedure are in agreement with 
those obtained using the worm algorithm and Ref. [45] . 
as shown in the main panel in Fig. Note that for the 
ground state calculation, lattice effects are negligible at 
the lowest densities studied. 


III. GROUND-STATE 
A. Homogeneous Systems 


1. Tonks-Girardeau limit 


In the Tonks-Girardeau limit, Tan’s contact allows for 
a simple analytic expression that unveils how to scale 
the momentum distribution function (one-particle corre¬ 
lations) to achieve data collapse for high momenta (small 
distances). As indicated in Eq. (10), C is linearly related 
to the interaction energy. In the language of Ref. |5| , 


C = 


NbPlvil) 

2-k 


( 11 ) 


where v{')) = p^')de(^)/d'^, and 6 ( 7 ) is a monotonically 
increasing function of 7 that saturates at 7 r ^/3 when 
7 —> 00 (for its definition, see Ref. [5]). For large val¬ 
ues of 7 , V has the following asymptotic form: 1 ^( 7 ) = 
4 p^e( 7)/(7 -I- 2). Thus, in the limit 7 —>• 00 , Tan’s con¬ 
tact reduces to C = 27rLp^/3. Hence, C = CjLp^ = 27r/3. 



k/p k/p 

FIG. 5. (Color online) Scaled momentnm distribution fnnc- 
tions m{k)/L in the Tonks-Girardeau limit as a fnnction of 
the scaled momentum k/p. (a) Systems with different number 
of particles and densities, (b) same as (a) but in log-log scale, 
(c) systems with the same number of particles and different 
densities, (d) systems with different nnmber of particles and 
the same density. Insets: (a) scaled average one-particle den¬ 
sity matrix, (c) and (d) same as main panels bnt in log-log 
scale. Thin dotted lines indicate the asymptotic behavior dis- 
cnssed in the text. Here, and throughout this work, p in the 
lattice is given in nnits of the inverse lattice spacing a. 


This means that if one plots m{k)/L vs k/p, no matter 
the number of particles and the system size, the curves 
must collapse for high values of k. From Eq. ([^, it then 
follows that the universal behavior of m{k) for high val¬ 
ues of k implies that, if one plots g{r)/p vs pr, the curves 
must collapse for small values of r. Results for m(fc) and 
g{r), obtained using the lattice approach discussed in the 
previous section, are reported in Fig. 

In Fig.ga), we plot m{k)/L vs k/p for systems with 
different number of particles and densities. Figure [^b) 
reports the same results but in a log-log scale. Both 
panels make apparent that, for intermediate and high 
momentum, all curves collapse onto a universal result. 
In the fc —^ 0 limit, m{k = 0)/L = mk=o oc y/Nb be¬ 
cause of the existence of quasi-long-range order | 2 ] and, 
as such, only curves with the same number of particles 
can collapse for low values of k. This can be better seen 
in Fig. j^c), where, for the same number of particles, all 
momentum distribution curves collapse for all momenta. 
Differences can only be seen at the highest momenta be¬ 
cause of finite-density effects in the lattice. On the other 
hand, as shown in Fig.j^d), plotting results for the same 
density and different number of particles, again produce 
curves that collapse at intermediate and high momenta 
but differ for the lowest values of k. For systems with 
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larger number of particles, the universal behavior can be 
seen to start at lower values of k. Consequently, irrespec¬ 
tively of the density in thermodynamically large system 
sizes, one expects the momentum distribution functions 
to become universal starting at infinitesimal values of k. 

In the inset in Fig.j^a), we plot results for g{r ) /p vs pr 
for systems with different number of particles and densi¬ 
ties. The plots can be seen collapse onto universal results, 
in which, at long distance, g{r) oc 1/^/r [5]. The latter 
behavior can be understood within the bosonization ap¬ 
proach [56], which predicts that g{r) oc at long 

distances, where K is the so called Tomonaga-Luttinger 
parameter. In the Tonks-Girardeau limit, K = 1. 


2. Finite interaction strength 

For finite interaction strengths, particles can pene¬ 
trate through one another and the systems tend to have 
lower kinetic energies (less particles in the high momen¬ 
tum tails) than in the Tonks-Girardeau limit. As a re¬ 
sult, Tan’s contact decreases with decreasing 7 . From 
Eq. ( |Tl] ), C can be written as C = L7^e'(7)p"‘/(27r). The 
scaled contact, C = CjLp'^ = 7^e'(7)/(27r), is now only 



FIG. 6 . (Color online) Scaled momentum distribution func¬ 
tions m(k)/L vs the scaled momentum k/p for finite values of 
7 , plotted in [(a) and (c)] a linear scale and in [(b) and (d)] a 
log-log scale. Results are reported for [(a) and (b)] 7 = 8.33 
and [(c) and (d)[ 7 = 1.0. Insets: g{r)/p vs pr for the same 
systems as in the main panels. In all plots, thin dotted lines 
indicate the asymptotic behavior discussed in the text. A 
linear fit of the long distance behavior of the curves in the 
insets indicate that the Luttinger parameter is A « 1.52 for 
7 = 8.33 and K ~ 3.63 for 7 = 1.0. Here, and throughout 
this work, p in the continuum is reported in an arbitrary unit. 


a function of 7 . From this result it follows that, so long 
as 7 is kept constant (no matter the system size and the 
number of particles in the system), the curves for m(k)/L 
vs k/p will be universal for large values of k. This is the 
same scaling discussed for the Tonks-Girardeau regime, 
in which 7 is constant (7 = 00 ). Similarly, the curves for 
g{r)jp vs pr must be universal at short distances. Results 
for m{k) and g{r), obtained using the worm algorithm^ 
are reported in Fig. 

In the main panels in Fig. we plot the momentum 
distribution for 7 = 8.3 [Fig. [^a) and [^b)] and for 
7 = 1.0 [Fig.j^c) and[^d)], in systems with different den¬ 
sities and average number of particles. Gonsistent with 
the previous discussion, the momentum distributions col¬ 
lapse for intermediate and high momentum. More im¬ 
portantly for experiments, with decreasing 7 , the onset 
of the k~^ behavior becomes less sharp (see also Fig. 
and Fig. 5 in Ref. [H]). This makes the experimental 
identification of the k~'^ tails increasingly difficult as one 
departs from the Tonks-Girardeau limit. 

As fc —>■ 0, the collapse obtained for intermediate and 
high values of k breaks down because of quasi-long-rage 
order. For finite values of 7 , more particles have low mo¬ 
menta than in the Tonks-Girardeau limit. In particular, 
within the bosonization approach, the zero momentum 
occupation is predicted to be m(k = 0 )/L a 
with A > 1 |5^. The insets in Figs.a) andj^c) depict 
g{r)/p vs pr for the same systems for which momentum 
distributions are shown in the main panels. The curves 
can be seen to collapse at short distances as expected. 
At long distances, g{r) oc as predicted by the 

bosonization approach, with K « 1.52 for 7 = 8.3 and 
K « 3.63 for 7 = 1.0. 


B. Trapped Systems 

1. Tonks-Girardeau limit 

In order to generalize the scaling discussed in 
Sec. |III ATj to harmonically trapped systems, we need to 
find the replacement for p and L in the presence of a trap. 
Let us start by discussing how to scale density profiles in 
a harmonic trap. Since the density profiles of impenetra¬ 
ble bosons are the same as those of the non-interacting 
spinless fermions to which they can be mapped | 2 |, we 
can focus on the noninteracting spinless fermions. The 
density profiles of the latter systems have been studied 
in detail in the past [57l - f5^ . For a homogeneous Fermi 
system, the density is a linear function of the Fermi mo¬ 
mentum kp, p = kp/p. The chemical potential p, in 
such systems is the Fermi energy Ep = h^k‘]p/{2m)^ so 
p = \J2mpl (Htt). 

In a harmonic trap, the local chemical potential 
changes according to the relation 

p{x) = po- V{x), 


( 12 ) 
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where is the chemical potential at the center of the 
trap. Within the local density approximation (LDA), the 
density at each position in the trap is solely determined 
by the local chemical potential, through the relation ob¬ 
tained for homogeneous systems. Hence 


n{x) 


1 

TT 



(13) 


where jl = 2m^jLo/h?. The integral of Eq. (131 over the 
entire space is the number of particles W- It provides the 
following relation between p, and Nt, p = 2 W/Qho- 
density in the center of the trap is then p = v^STYb/TraHO j 
and the position in the trap at which the density vanishes 
is L = ^/2N},aB_o- Using these quantities, the density 
profile in the trap can be written as 


n{x) = p 


1 - 



(14) 


By defining the scaled density h{x) = n{x)/p and the 
scaled position x = xjL, one obtains 

h{x) = (1 — i^)5. (15) 



FIG. 7. (Color online) [(a) and (c)] Universal scaling of den¬ 
sity and [(b) and (d)] momentum distribution functions in the 
Tonks-Girardeau limit. We report results from lattice calcu¬ 
lations for [(a) and (b)] up to 25 bosons, and from the worm 
algorithm for [(c) and (d)] up to 100 bosons. In (a) and (c), 
we also report the LDA prediction for the density profiles. In 
(b) and (d), thin dotted line depict k~'^ behavior. 


Since the momentum distribution function at high mo¬ 
menta is determined by the average short-distance one- 
particle correlations, from the LDA and Eq. (15) one can 
advance that, using the scaling relations for homogeneous 
systems with p ^ p and L ^ there will be data col¬ 
lapse for the momentum distribution function at high 
momenta. 


In Fig. we report results for density and momen¬ 
tum distribution functions in the Tonks-Girardeau limit 
obtained using the lattice approach and the worm algo¬ 
rithm. Given the constraint of very low site occupancies 
in the lattice, our lattice calculations can only be done 
for significantly smaller numbers of particles than within 
the worm algorithm. The normalized density profiles in 
Figs. [3a) and [3 c) exhibit the expected collapse, even for 
as few as nine bosons [Fig. [^a)]. Data collapse can also 
be seen in the momentum distribution functions, except 
at very low momenta. The k~'^ high momentum tails are 
apparent both in the lattice and worm algorithm results, 
though they are significantly better seen in the lattice cal¬ 
culation where much smaller values of m{k), for higher 
values of k, can be computed. 


We note that, using the relations we established in 
Sec. |ll| between quantities in the continuum and in the 
lattice, one finds that the scaling properties of trapped 
systems in the continuum discussed so far are consistent 
with those discussed for lattice systems in Ref. [HI |33] . 


2. Finite interaction strength 


For finite interaction strengths, the relation between 
the chemical potential and the density is not as simple as 
in the Tonks-Girardeau limit, where p, = p^/{2m). 

Instead, p = p'^ f (j)/{2m), where /(y) = 3 e( 7 )— 7 e'( 7 ) 
[5]. This means that, using Eq. (12) and within LDA, 


n\x)f[^{x)]=n^{0)f[j{0)]-^. 


(16) 


^HO 


The position in the trap at which the density vanishes 
is Ull = ^/[7(0)]n(0)a^o. Defining, as for the Tonks- 
Girardeau limit, n{x) = n{x)/n{0) and x = x/Lll, and 
noticing that 7 ( 0 ;) = mg/[h‘^n{x)] = 'y{0)/n{x), we ob¬ 
tain 


n^{x)f 


7(0) ' 

n{x) 


f[j{0)]{l-x^). 


(17) 


Equation ( jlTj ) makes apparent that the scaled density 
profiles of systems with the same value of 7 in the center 
of the trap must be identical. As in the Tonks-Girardeau 
limit, this result and the LDA allow one to advance a 
scaling collapse in the momentum distribution of trapped 
systems at high momenta [28j . 

However, the scaling above requires the knowledge of 
n( 0 ) and All, which can only be obtained numerically 
(or measured experimentally). We can predict another 
universal way to scale densities and positions by noting 















FIG. 8 . (Color online) [(a) and (c)] Universal behavior of 
scaled density profiles and [(b) and (d)] momentum dis¬ 
tribution functions for harmonically trapped systems with 
7 ( 0 ) « 9.3 in (a) and (b), and 7 ( 0 ) « 0.5 in (c) and (d). 
In (b) and (d), thin dotted lines depict k~'^ behavior. 


that, from Eq. ( [T7| ), it follows that 

h(f) =/io[l - ^^ 7 ( 0 )], (18) 

where hQ(x,y) is an unknown function that can, in 
principle, be computed numerically (from now on 
we will call dummy functions of this kind as ha, 
with a = 0,1,...). Since = 2 n(x)dx = 
2n^(0)a^oVfh(0)] fo holl - x^,j(0)]dx, we see that 
n( 0 ) = hi[ 7 ( 0 )]'\/]Vh/aHO = ^ 2 [ 7 ( 0 )]p, and p is the den¬ 
sity in the center of the trap in the Tonks-Girardeau limit. 
From this result for n(0) it follows that Lll = ^ 3 [ 7 ( 0 )]E, 
where L is the point at which the density vanishes in 
the Tonks-Girardeau limit. Hence, after fixing 7 ( 0 ), one 
can use exactly the same scaling for finite interaction 
strengths as the one used in the Tonks-Girardeau limit. 
We should stress that this is unique to bosons in the 
continuum. In the Bose-Hubbard model, one needs two 
parameters (the so-called characteristic density and the 
on-site interaction strength) to be kept fixed in order to 
be able to scale density profiles m- 

In Fig. 1^ we show the scaled density [Figs. |^a) and 
| 8 ][c)] and momentum [Figs, [^b) and |^d)] profiles for 
systems with different numbers of particles, trapping fre¬ 
quencies, and ID scattering lengths. Results are pre¬ 
sented for 7 ( 0 ) « 9.3 [Figs.j^a) and[^b)] and 7 ( 0 ) « 2.9 
[Figs.§c) and l^d)]. The plots exhibit an excellent data 
collapse in the density profiles, and in the momentum dis¬ 
tribution functions at intermediate and high momentum, 
as predicted by the analysis above. The high momen¬ 
tum k~'^ tails are barely apparent in Figs.|^b) and[^d). 


By comparing the results for 7 ( 0 ) « 9.3 to those in the 
Tonks-Girardeau limit in Fig. one can see that, even 
for that large value of 7 ( 0 ), the range and sharpness of 
the k~‘^ tails have decreased. This worsens as 7 ( 0 ) is fur¬ 
ther decreased and, as in the homogeneous case, poses an 
increasing experimental challenge for detecting the 
tails as one departs from the Tonks-Girardeau regime. 


IV. FINITE TEMPERATURE 
A. Homogeneous Systems 

1. Tonks-Girardeau limit 

At finite temperature, the long-distance behavior of 
one-particle correlation functions is qualitatively differ¬ 
ent from the one in the ground state. This because 
quasi-long-range order is destroyed by thermal fluctua¬ 
tions, i.e., the correlation functions decay exponentially 
with increasing distance [2]. As a result, the population 
of all momentum modes m{k)/L becomes intensive. This 
means that, if the density and the temperature are kept 
constant, the momentum distribution function of systems 
with different number of particles should collapse (pro¬ 
vided the system sizes are much larger than the corre¬ 
lation length). The curves with the same value of p in 
Figs, [^a) and[^b) show that this is indeed the case. 

The question we would like to address here is how to 
scale the momentum profiles when the density and the 
temperature are changed in these systems. For that, as 
we did in Sec. ImATI we look into Tan’s contact in the 
limit 7 —)■ 00 . As follows from the derivation in Ap- 



FIG. 9. (Golor online) Finite temperature results for the 
scaled momentum distribution function of the homogeneous 
Tonks-Girardeau gas in systems with different average num¬ 
ber of particles and densities. The scaled temperature in the 
plots is: (a) T = 2.51 and (b) f = 15.06. The black solid 
line in both panels was obtained using the lattice approach, 
while all other curves were obtained using the worm algo¬ 
rithm. Thin dotted lines depict behavior. 
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pendixj^ the Tan’s contact takes the form 

\ r / M 


_ {2mkBT)'^L 
^ ~ h^2TT 


kuT 


f3/2 


kuT 


(19) 


where fi,{ ) is the Fermi-Dirac function (see Appendix [C| 
for its definition). From Eq. (C7), one can see that 
^/(fcsT) is a function of 2mkBT/p^). Thus, one 
can define the scaled (dimensionless) temperature T = 
2mkBT/{h^p^). The scaled contact C = C/{Lp‘^) then 
has the form 




/3/2 




( 20 ) 


which is only a function of T. This means that, once 
T is fixed, m{k)/L vs k/p is universal. The results in 
Figs. ia) and [^b) , in which we plot systems with dif¬ 
ferent densities and temperatures but the same values 
of T, show that this is indeed the correct scaling. The 
momentum tails are clearly visible at high momenta. 
Two important features of those tails are apparent when 
comparing Figs, [^a) and Sb). The first one is that the 
weight of the tails (the contact) increases with increasing 
T, as discussed in Ref. [45]. The second one, probably 
more important for experiments, is that with increasing 
T the values of m{k) at which the tails start to develop 
become smaller. This means that the k~‘^ momentum 
tails are more difficult to observe experimentally as the 
temperature increases, despite the fact that Tan’s con¬ 
tact becomes larger. 


2. Finite interaction strength 

For systems with a finite value of 7 and finite temper¬ 
atures one could advance, based on the results for the 



FIG. 10. (Color online) Finite temperature results for the 
scaled momentum distribution function of the homogeneous 
Lieb-Liniger gas with different average number of particles 
and densities at finite 7 = 8.33. The scaled temperature in 
the plots is: (a) T = 2.51 and (b) T = 15.06. All results were 
obtained using the worm algorithm. Thin dotted lines depict 
k~'^ behavior. 


ground state (in which scaling collapse was found so long 
as 7 was kept fixed) and on the results for finite tem¬ 
perature in the Tonks-Girardeau limit (in which scaling 
collapse was found so long as T was kept fixed), that scal¬ 
ing collapse requires us to keep 7 and T constant. The 
results reported in Fig. W for a finite value of 7 = 8.33, 
which parallel the results in the Tonks-Girardeau limit 
reported in Fig.j^ show that this is indeed the case. The 
behavior of the high momentum tails with increasing T 
in Fig. is qualitatively similar to that discussed in the 
Tonks-Girardeau limit. 


While a close expression for Tan’s contact for finite 
values of 7 and finite temperatures is not available, a 
calculation up to order I /7 (see Appendix [Cj) reveals that 
C = Cj (Tp^) can be written as [6TVI63] 


C = 


f2 

27r 


/3/2 X 


•^ 1/2 + ^ (2/1/2 + /-I/2/3/2) 


, ( 21 ) 


where by fi, it is meant /^{p/ksT). This further sup¬ 
ports the correctness of the scaling used in Fig. [T^ 



FIG. 11. (Golor online) Finite temperature results for [(a) 
and (c)] the scaled density profiles and [(b) and (d)] momen¬ 
tum distribution functions of harmonically trapped Tonks- 
Girardeau gases with [(a) and (b)] T{0) ~ 2.63 and [(c) and 
(d)[ T'(O) « 37.48. In (a) and (c), black solid lines show the 
LDA predictions. In (b) and (d), the black solid lines show 
results of the lattice calculations. All other curves were ob¬ 
tained using the worm algorithm. In (b) and (d), thin dotted 
lines depict k~'^ behavior. 
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B. Trapped Systems 

1. Tanks-Girardeau limit 


As for the ground-state case, let us discuss how to 
scale density profiles in the Tonks-Girardeau limit at fi¬ 


nite temperature. Dividing Eq. (12) by we obtain 


^ _ 

keT keT 2kBT ’ 


/r(0) muj'^x^ 


( 22 ) 


Using the LDA and the fact that, as discussed in the 
previous section, in homogeneous systems ^/(ksT) = 
hi{2mkBT/h?p'^), we can rewrite Eq. (22) as 


hi 


f(0)n2(0) 

n^(x) 


= h4[fm - 


a%on0)n^i0y 


(23) 


where T(0) = 2mkBT/[h‘^n^{0)]. Now, let us find the 
position Lpx in the trap at which the chemical poten¬ 
tial vanishes. At zero temperature, the position at which 
the chemical potential vanishes is also the position at 
which the density vanishes. However, this is not the 
case anymore at finite temperature. For p{x) = 0, 
/i 4 [ 2 mA:BT//i^n^(a;)] = 0 and, using Eq. (23), Lft = 

\jT{0)h4 [r(0)]n(0)a|^Q. Defining x = x/L^t, Eq. (23) 
can be rewritten as 


hi 


f(0)n(0f 

n{x)‘^ 


= hi[f{Q)][l-i^). 


(24) 


Defining fi{x) = n(x)/n(0), we see that Eq. (24) implies 
that 


fi{x) = h^ll - i^,r( 0 )]. 


(25) 


At this point, we note that Eq. (25) has exactly the same 
form as Eq. ( [I^ , with ho ho and 7 ( 0 ) —>■ T’(O). Using 
that Nb = 2 n{x)dx, we get that n( 0 ) = hQ\f{Q)]p 
and Lft = h'!\T{Q)]L. This means that, once 'r(O) has 
been fixed, a scaling as the one used for the ground state 
case will produce universal density profiles at finite tem¬ 
perature. In other words, for a fixed T(0), p and L deter¬ 
mine the average density in the system and an effective 
system size. Hence, given the fact that LDA is expected 
to work for sufficiently large system sizes at any finite 
temperature (because of finite correlation lengths), plots 
of m{k)/L vs k/p are expected to exhibit data collapse. 

In Figs. [II|(a) and ©c), we show plots of n{x)/p vs 
x/L for systems with different average number of parti¬ 
cles, densities in the center of the trap, and temper atur es, 
but all with the same value of 'r(O) 
and f (0) « 37.48 [Fig. 
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2.63 [Fig. 11 a)] 


c)]. All the scaled density pro¬ 
files exhibit the predicted data collapse and are in excel¬ 
lent agreement with the LDA prediction. In Figs. ©b) 
and ©d), we show the corresponding scaled momen- 
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FIG. 12. (Color online) Finite temperature results for [(a) and 
(c)[ the scaled density profiles and [(b) and (d)[ momentum 
distribution functions of harmonically trapped Lieb-Liniger 
gases with [(a) and (b)[ 7 ( 0 ) « 9.3, and T'(O) « 1.56 and [(c) 
and (d)] T’(O) « 21.38. All curves were obtained using the 
worm algorithm. In (b) and (d), thin dotted lines depict k~‘^. 


turn distribution functions, which also exhibit an almost 
perfect data collapse. The results for the momentum 
distribution functions exhibit an excellent agreement be¬ 
tween the worm algorithm calculations and the lattice 
approach results. The latter allow us to reach smaller 
values of m{k)/L (larger values of k/p) so that the k~^ 
momentum tails can be better resolved. 


2. Finite interaction strength 


The scaling of trapped systems with finite interaction 
strengths at finite temperature is the most challenging 
one. From the ground state analysis in the presence of 
a trap we know that 7 ( 0 ) needs to be kept fixed. From 
the finite temperature analysis of homogeneous systems 
we learned that 7 and T need to be fixed, while the finite 
temperature analysis of the trapped Tonks-Girardeau 
gas revealed that T’(O) needs to be fixed. Putting all 
this together we can advance that, for finite interaction 
strength, trapped systems at finite temperature will only 
exhibit scaling collapse if 7 ( 0 ) and T(0) are fixed. We can 
also advance, from the analysis of trapped systems with 
finite interaction strength in the ground state and in the 
Tonks-Girardeau limit at finite temperature that using p 
and L will allow one to achieve the expected collapse. 

In Fig. we show results for the scaled momentum 
distribution function of trapped s_ystems with the same 
value of 7 ( 0 ) and two values of 'r(O). They all exhibit 
a nearly perfect collapse. The features observed in the 
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momentum tails of the homogeneous case are also appar¬ 
ent here. Namely, Tan’s contact increases with increasing 
temperature but at the same time the tails develop 
starting at higher values of k and smaller values of m{k). 

V. CONCLUSIONS 

We carried out an unbiased study of one-dimensional 
trapped bosons in the ground state and at finite tem¬ 
perature using the worm algorithm, and the Bose-Fermi 
mapping in the presence of a lattice at low fillings. When¬ 
ever possible, we compared our results to the predictions 
of the Bethe ansatz obtaining an excellent agreement. 

Our study focused in the behavior of the density and 
momentum distribution functions, as well as one-particle 
correlations. We discussed in detail how to scale density 
and momentum profiles to observe universal behavior, 
which was demonstrated by our numerical results. For 
trapped systems in the ground state, we showed that fix¬ 
ing 7 in the center of the trap is all that is required to 
obtain a universal scaling of density and momentum pro¬ 
files. This is to be contrasted with the Bose-Hubbard 
model, for which one needs to fix the so-called character¬ 
istic density and the on-site interaction strength . At 
finite temperature in systems in the continuum, in ad¬ 
dition to the condition for the ground state, one needs 


to fix the scaled temperature in the center of the trap. 
Those two parameters fully characterize the densities and 
momentum distribution functions. 

In our study of the momentum distribution function, 
we payed special attention to the k~'^ asymptotic be¬ 
havior at high momenta. As shown by our numeri¬ 
cal calculations, those momentum tails become increas¬ 
ingly challenging to resolve as the temperature increases 
and the interactions become weaker, making the Tonks- 
Girardeau limit at low temperatures the ideal regime for 
their experimental observation. 
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Appendix A: Pair-product action 


The pair-product action assumes that interactions be¬ 
tween different pairs of particles are uncorrelated. Thus, 
the action from the interaction breaks into sum of pair 
actions 


U2{R,R'\t) =^U2{rij,r[j]T). (Al) 

i<j 

The pair action for particles with repulsive contact in¬ 
teractions can be obtained exactly from the relative two- 
body propagator ps, U 2 = — Inp^. The relative two- 
body propagator is the two-body propagator divided by 
the free propagator ps = G 2 IGQ. For the Lieb-Liniger 
Hamiltonian, written as 


N N 

H = g ^ 5{x,-xj), 

i-l i<i = l 


(A2) 


the relative propagator can be written as 


\ / — X eric 

>1 + \y\ +9x' 

V A' 

•\/4A'r 


X exp 


^ {x-yf + ^{\x\ + \y\)+'^ 


4A't 


4A' 


(A3) 


where A' = 2A and erfc[x] = dtexp{—t‘^) is the 

complementary error function (see Refs. |64l I65j for the 
analytic derivation of this result). 


Appendix B: Energy estimators 


The thermodynamic estimator for the total energy 
is computed from Et = —Z~^dZ/dP, while the 
one for the kinetic energy is computed from Kt = 
—m/3~^Z~^dZ/dm [5T1|53]. One obtains that 


Ej^ = 


Kx = 


Vt = 



{ri,-r:e^)^ 

4At^ 

(jIq, -TIq, j 

4Ar2 

XdU^\ 

X dX X 


/a., 

A^\ 
T dX / 


(Bl) 


Here, i and a are indices for particle and discrete imag¬ 
inary time, respectively. and its derivatives can be 
computed using the pair-product action discussed in Ap¬ 
pendix 
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Appendix C: Yang-Yang Thermodynamics 


The Yang-Yang equation for the dressed energy is 
=TiZ- - M 


2 m 
keT 
27r 


dq 


2 c 


+ {k — qY 


ln[l 


<i) 

■ g 


(Cl) 


e(fc) is used to compute thermodynamic quantities. The 
thermodynamic potential density is given by the expres¬ 
sion 


n{n,c,T) = - 


keT 

2tt 


dq ln[l 


- g 


(C2) 


All other thermodynamic quantities of interest here can 
be derived from n(/i, c, T) through the following relations 

e = il + fj.p + Ts, p = -ll, p = - — , s = - — , (C3) 

where e is the energy density, s is the entropy density, 
and p is the pressure. 


For very strong contact interactions, i.e., c ^ 1, one 
can Taylor expand e(fc) in terms of 1/c 

The dressed energy e{k) can then be obtained iteratively. 
To zeroth order: = h^k^/{2m) — p. This allows 

one to compute Q. to lowest order 

^ , (C5) 


where 


fy{<y) = 


r(^) 


g(a:-a) _|_ l' 


(C6) 


is the Fermi-Dirac function. Using Eq. (C5), and the 


relations (C3), one obtains the following expressions for 


the density and the energy of the system 

( 0 ) j: ( M 

" = \k^ 

JO) ^ ( P 


(C7a) 

(C7b) 


Substituting fl —)• in Eq. (C4), one obtains 


= k^ — p^^\ where 
, s/^{kBT)'^/'^ j. 

P^ ’ = P+ . f3/2 , , ^ 

hy/TTC 


(C8) 




Now, substituting p -J p^^'> in Eq. (C5), and expanding 
to first order in 1/c, we get 

n<‘> = X (1 + y:^A, 

hj2TT \ chYn 

(C9) 

where, from now on, by fi, it is me ant fuip/ksT). Fi¬ 
nally, substituting fl -J in Eq. (C4|, we get that 


(CIO) 


Y‘^\k) = k^ — p^'^\ where 

(2) , V^{kBT)^^‘^ ^ 

- M + r- ^3/2 

AiyTTC 

/ y/2.mkBT 

" (' + 

This expression allows us to obtain by substituting 
p -J p^‘^'> in Eq. (C5). Expanding the resulting equation 


up to second order in 1/c, one obtains 
0 ( 2 ) v^(fcBr)3/2 |- ^2mkBT 


( \/2mkBT 

\ ch\pk 


(^fi/2 + 2 / 3 / 2 /- 1 / 2 ^ 


• (Cll) 


Using Eq. (|C11|), we can calculate Tan’s contact 

(C12) 


through Tan’s sweep relation 

2m(YL / dp 


C = -- 


S 2 


dc 




The result obtained to order 1/c is: 
_ {2mkBT)‘^L 


(C13) 


v^(fcBr)i /2 2 . . X 

- nX /^ c —/-1/2/3/2) 


These are nothing but the density and the energy density 
of a system of noninteracting spinless fermions. 












































